Escolhendo a largura de banda

Vamos agora explorar, na prática, como diferentes larguras de banda na Estimativa de Densidade de Kernel afetam a suavização quando nossos dados estão espacialmente clusterizados. No R, podemos testar vários valores com a função density(). Essa função possui alguns parâmetros do Kernel como:

• Largura de banda (sigma): Define o alcance de suavização do kernel.
• Tipo de kernel (kernel): Permite escolher o formato do kernel.
• Correção de bordas (diggle): Ajusta a densidade para considerar bordas da janela de observação.
• Vetor de pesos (weights): Define pesos diferentes para os pontos.

Neste caso, vamos alterar o valor de sigma na função density() utilizando o objeto agregado que foi criado para os exemplos anteriores. Acompanhe o código abaixo:

kernel_cluster1 <- density(agregado, sigma = 0.03)
kernel_cluster2 <- density(agregado, sigma = 0.08)
kernel_cluster3 <- density(agregado, sigma = 0.15)

# vamos apresentar agora os gráficos em ggplot
library(ggplot2)
# vamos carregar a biblioteca patchwork para
# plotar os gráficos juntos
library(patchwork)


g_den1 <- 
  ggplot(as_tibble(kernel_cluster1), aes(x, y)) +
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(agregado),
             color = "white") + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Largura de banda = 0,03") +
  guides(fill = "none") + 
  theme_void()

g_den2 <- 
  ggplot(as_tibble(kernel_cluster2), aes(x, y)) +
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(agregado),
             color = "white") +
  scale_fill_viridis_c(option = "B") +
  labs(title = "Largura de banda = 0,08") +
  guides(fill = "none") + 
  theme_void()

g_den3 <-
  ggplot(as_tibble(kernel_cluster3), aes(x, y)) +
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(agregado),
             color = "white") + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Largura de banda = 0,15") +
  guides(fill = "none") +
  theme_void()

(g_den1 | g_den2 | g_den3)

A Figura 10 mostra três mapas de calor gerados pelo método de estimativa de densidade de Kernel, aplicados a diferentes larguras de banda ($\tau$) representadas pelos valores 0,03, 0,08 e 0,15. Cada painel corresponde a uma largura de banda específica:

• Largura de banda = 0,03: O mapa apresenta alta sensibilidade a variações locais, com áreas de concentração muito definidas e detalhadas. A suavização é mínima, evidenciando pequenos agrupamentos ou hotspots de alta densidade.

• Largura de banda = 0,08: O aumento da largura de banda reduz a granularidade da análise, resultando em uma distribuição mais suavizada. Os padrões locais começam a se mesclar, destacando áreas de concentração maiores, mas com menos detalhes.

• Largura de banda = 0,15: Com a maior largura de banda, a suavização é intensa, e os agrupamentos menores se fundem em grandes áreas de alta densidade. Este nível de suavização destaca padrões globais, mas pode ocultar variações locais.

Resumindo, a Figura 10 demonstra como a largura de banda no método Kernel afeta a suavização espacial. Larguras menores enfatizam detalhes locais, enquanto larguras maiores evidenciam padrões globais. A escolha da largura de banda deve equilibrar a necessidade de “capturar” detalhes locais sem perder a “visão geral” da distribuição dos eventos.

O pacote ggplot2 também possui funções próprias que permitem a visualização de estimativas de densidade, incluindo o uso de Kernel espacial. Uma dessas funções é a stat_density2d(), que ajusta o Kernel diretamente na visualização e adiciona automaticamente as linhas de contorno, facilitando a interpretação das áreas de maior ou menor densidade. A largura de banda mais adequada para a estimação do Kernel é determinada automaticamente, através do pacote MASS. Abaixo, vamos comparar a largura de banda “ótima” com uma largura de banda arbitrária, escolhida sem muitos critérios. Acompanhe abaixo:

g_contour <- 
  ggplot(as_tibble(agregado), aes(x, y)) +
  stat_density2d_filled(h = c(1, 1)) + 
  stat_density2d(h = c(1, 1), n = 200, contour_var = "count") + 
  geom_point(color = "white") +
  labs(title = "Largura de banda arbitrária (0,1)") + 
  guides(fill = "none") +
  theme_void() 

# vamos inserir linhas de contorno também

opt_bw <- c(MASS::bandwidth.nrd(agregado$x), MASS::bandwidth.nrd(agregado$y))

g_contour_opt <- 
  ggplot(as_tibble(agregado), aes(x, y)) + 
  stat_density2d_filled(h = opt_bw) + 
  stat_density2d(h = opt_bw, n = 200, contour_var = "count") + 
  geom_point(color = "white") +
  labs(title = "Largura de banda ótima") +
  guides(fill = "none") + 
  theme_void()

(g_contour | g_contour_opt)

Na Figura 11 podemos observar que a largura de banda otimizada reflete a distribuição real dos dados com maior precisão, enquanto a largura arbitrária pode gerar resultados mais espúrios.

Mesmo com uma largura de banda bem escolhida, suas estimativas podem ser subestimadas nas periferias do mapa. Por isso, no item seguinte veremos correções de borda, muito útil quando a área de estudo é limitada (por exemplo, um bairro cercado por barreiras geográficas).

Estimativa de Kernel com correção por bordas

Como citado anteriormente, o objetivo da correção por bordas na estimativa de densidade de Kernel é lidar com o viés introduzido nas extremidades da região de estudo. Ou seja, em áreas próximas às bordas, o interpolador Kernel pode saltar para fora da região de interesse, subestimando a densidade real, pois há menos pontos dentro do raio de influência do Kernel comparado ao que ocorre no interior da região.

A correção por bordas lida com esse problema ao considerar apenas a contribuição efetiva do Kernel dentro da região de estudo. Isso é feito normalizando o volume do Kernel sobre a área que está realmente dentro da região.

Para calcular essa estimativa, primeiramente, calculamos o volume sob o Kernel que está efetivamente dentro da região de estudo, representado por:

$$\delta_\tau(s) = \int_R \frac{1}{\tau^2} k\left(\frac{s - u}{\tau}\right) du$$

Onde:

• $\delta_\tau(s)$ é o volume sob o Kernel ajustado pela região de interesse.

• $k(\cdot)$ é a função Kernel escolhida.

• $s - u$ é a diferença entre a localização de interesse ($s$) e os pontos dentro da região ($u$).

• $\tau$ é o parâmetro de largura de banda que controla a suavização.

• $R$ é a região de estudo.

Após calcular essa correção, o estimador de densidade ajustado é dado por:

$$\hat{\lambda}(s) = \frac{1}{\delta_\tau(s)} \sum_{i=1}^n \frac{1}{\tau^2} k\left(\frac{s - s_i}{\tau}\right)$$

Onde:

• $\hat{\lambda}(s)$ é a estimativa da densidade corrigida.

• $\delta_\tau(s)$ é o fator de correção do Kernel calculado anteriormente.

• $n$ é o número total de pontos observados.

• $k(\cdot)$ é a função Kernel.

• $s - s_i$ é a diferença entre a localização de interesse ($s$) e os pontos observados ($s_i$).

• $\tau$ é o parâmetro de largura de banda.

Quando falamos em correção por bordas e diversas fórmulas, pode parecer complicado, mas não se preocupe! A correção por bordas é feita acrescentando o parâmetro diggle = TRUE na função density(), que já vimos. Agora, vamos comparar a função de Kernel utilizando o padrão de pontos agregado e a mesma largura de banda com e sem a correção por bordas. Veja o código abaixo:

kernel_cluster1 <- density(agregado, sigma = 0.08)
kernel_cluster2 <- density(agregado, sigma = 0.08, diggle = TRUE)

# vamos apresentar agora os gráficos em ggplot

g_den1 <- 
  ggplot(as_tibble(kernel_cluster1), aes(x, y)) +
  geom_tile(aes(fill = value)) +
  geom_point(data = as_tibble(agregado),
             color = "white") + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Sem correção por bordas",
       subtitle = "Largura de banda = 0,08") +
  guides(fill = "none") + 
  theme_void()

g_den2 <- 
  ggplot(as_tibble(kernel_cluster2), aes(x, y)) +
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(agregado),
             color = "white") + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Com correção por bordas", 
       subtitle = "Largura de banda = 0,08") +
  guides(fill = "none") + 
  theme_void()

(g_den1 | g_den2)

Observe na Figura 12 que a correção por bordas aumentou a intensidade nos pontos na borda inferior e lateral direita. A correção por bordas tende a ajustar a densidade para compensar o viés nas bordas da área de estudo, resultando em estimativas mais precisas para regiões próximas aos limites.

Viu como foi fácil? Em seguida, vamos explorar as diferenças entre tipos de interpolador Kernel e como isso pode impactar a análise.

Diferenças entre as funções de Kernel

A função Kernel é um componente essencial na estimativa de densidade de Kernel, pois define como os pontos vizinhos influenciam a estimativa de densidade em um ponto específico. Mas, diferentes funções de Kernel podem resultar em superfícies de densidade distintas, mesmo com a mesma largura de banda. Muitos softwares (como o QGIS, por exemplo), padronizam o uso de uma só função, apesar que disponibilizar outras opções. Neste curso, veremos três tipos de Kernel: Gaussiano, Quártico e o Disc. É importante entender suas diferenças e o impacto na estimativa de densidade, mesmo com a mesma largura de banda.

Para visualizarmos as diferenças entre as funções de Kernel, vamos aplicar três tipos diferentes de Kernel ao mesmo conjunto de dados (o padrão agregado) e com a mesma largura de banda. Acompanhe o código abaixo e, em seguida, vamos discutir os resultados.

kernel_cluster1 <- density(agregado, kernel = "gaussian", sigma = 0.08)
kernel_cluster2 <- density(agregado, kernel = "quartic", sigma = 0.08)
kernel_cluster3 <- density(agregado, kernel = "disc", sigma = 0.08)

# vamos apresentar agora os gráficos em ggplot

g_den1 <-
  ggplot(as_tibble(kernel_cluster1), aes(x, y)) +
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(agregado),
             color = "white") + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Kernel Gaussiano", 
       subtitle = "Largura de banda = 0,08") +
  guides(fill = "none") + 
  theme_void()

g_den2 <- 
  ggplot(as_tibble(kernel_cluster2), aes(x, y)) +
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(agregado),
             color = "white") + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Kernal Quártico", 
       subtitle = "Largura de banda = 0,08") +
  guides(fill = "none") + 
  theme_void()

g_den3 <- 
  ggplot(as_tibble(kernel_cluster3), aes(x, y)) +
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(agregado),
             color = "white") +
  scale_fill_viridis_c(option = "B") +
  labs(title = "Kernel Disc",
       subtitle = "Largura de banda = 0,08") +
  guides(fill = "none") + 
  theme_void()

(g_den1 | g_den2 | g_den3)

Observando a Figura 13, podemos notar diferenças importantes entre os tipos de Kernel. Vamos analisar cada um deles:

• O Gaussiano utiliza uma distribuição normal, o que significa que os pontos próximos ao centro do Kernel recebem pesos maiores, mas o decrescimento do peso ocorre de forma contínua e suave à medida que a distância aumenta. As transições de cores são suaves e graduais, refletindo a suavização contínua típica do Kernel Gaussiano. Isso resulta em áreas de maior densidade (em laranja/vermelho) são bem definidas, mas as bordas das áreas têm transições gradativas, o que cria uma aparência mais difusa.

• O Quártico (também conhecida como Epanechnikov) dá mais peso a pontos próximos do que a pontos distantes, mas o decrescimento do peso é gradual e termina abruptamente em um limite definido (como uma “bolha”). As áreas de alta densidade permanecem bem definidas, mas as transições para áreas de menor densidade são ligeiramente menos difusas que no Kernel Gaussiano. Cria uma suavização intermediária, que destaca os agrupamentos sem perder muita precisão local.

• O Disc dá peso igual a todos os pontos dentro de um raio fixo e ignora pontos fora desse raio. Isso cria áreas de densidade com contornos abruptos e sem transições suaves. As áreas de alta densidade aparecem como “manchas” com bordas bem definidas, sem transições graduais para áreas de menor densidade. Esse Kernel é ideal para destacar claramente os agrupamentos locais, mas sacrifica a suavidade dos resultados.

Mas, no dia a dia da vigilância em saúde, qual função Kernel devemos escolher? Depende do objetivo! Aqui estão algumas orientações:

• Gaussiano: se o objetivo for explorar padrões gerais de distribuição dos pontos, o Kernel Gaussiano pode ser mais adequado, pois destaca padrões gerais em grandes áreas.

• Quártico: se o objetivo for identificar agrupamentos locais, com dados clusterizados, o Kernel Quártico pode ser mais apropriado, pois mantém detalhes locais e balanceia a suavização nas áreas mais distantes.

• Disc: se o objetivo for detectar áreas “hotspots”, com agrupamentos nítidos, o Kernel Disc pode ser mais apropriado, pois destaca claramente os pontos de alta densidade.

Essas escolhas permitem ajustar a suavização aos objetivos específicos da análise espacial. Mas, a análise utilizando estimadores de Kernel não se limita apenas a suavizar a distribuição dos eventos. Também é possível incorporar atributos e calcular razões entre diferentes eventos, o que pode enriquecer ainda mais a análise espacial. Vamos ver como isso funciona na prática nos próximos itens.

Kernel por atributo (kernel ponderado)

Na rotina da vigilância, nem sempre basta saber onde ocorrem os casos; é crucial relacioná-los a um denominador: por exemplo, população exposta, leitos disponíveis, número de criadouros, entre outros. A Estimativa de Densidade de Kernel resolve isso facilmente: basta atribuir um peso a cada ponto, representando o valor do atributo de interesse. Assim, a estimativa de densidade ponderada é calculada como uma média ponderada dos pontos vizinhos, levando em conta o valor do atributo definido.

Neste sentido, o método de suavização utilizando o Kernel permite incorporar uma covariável (atributo) para refinar a estimativa. Por exemplo, é possível estimar a população por unidade de área ou calcular a razão entre dois Kernels, gerando uma estimativa suavizada de eventos por população.

A fórmula para a estimativa por atributo ($y_i$) é:

$$\hat{\lambda}(s) = \sum_{i=1}^{n} \frac{1}{\tau^2} k \left( \frac{s - s'_j}{\tau} \right) y_i$$

Onde:

• $\hat{\lambda}(s)$ é a estimativa do atributo para unidade de área.
• $\tau$ é a largura de banda.
• $y_i$ é o valor do atributo associado a cada ponto $i$.
• $k$ é a Função Kernel, que define como os valores dos pontos vizinhos contribuem para a estimativa.

Este tipo de análise permite representar espacialmente tanto a distribuição do atributo (ex.: população) quanto a intensidade de eventos proporcional ao atributo, fornecendo uma visão mais detalhada e ajustada à realidade.

Quadro 3: Exemplo de aplicações de Kernel por atributo.

Atributo	Como	Pergunta	Aplicação
População do setor censitário (habitantes)	Atribuindo o número de habitantes ao centroide ou ao centro populacional da área	Onde se concentram mais pessoas?	Planejar vacinação, larvicida
Eventos / População	Calculando a razão de eventos (ex.: casos de doença) em relação à densidade populacional.	Qual é o risco por 1.000 hab.?	Comparar bairros ajustando pelo tamanho populacional

Mas, como isso funciona na prática? O próximo passo agora é combinar duas dessas superfícies: dividir o kernel de casos pelo kernel de população (ou por outro denominador relevante) para obter uma razão de kernels, isto é, um mapa suavizado de risco ou taxa. No próximo item vamos ver, na prática, como gerar essas duas camadas no R e calcular a razão de forma rápida, usando as mesmas funções que acabamos de comentar.

Razão de Kernel

Até aqui, vimos como gerar uma superfície suavizada para qualquer conjunto de pontos. Na rotina da vigilância, porém, o que mais interessa em determinados momentos é a taxa. Por exemplo, quantos eventos ocorrem para cada habitante (ou para cada leito, criadouro, domicílio visitado, entre outros).

A razão de Kernel faz exatamente isso. É uma técnica utilizada para calcular uma taxa suavizada a partir da divisão de duas estimativas de Kernel. Esse método é especialmente útil em análises espaciais que relacionam eventos com outros eventos (por exemplo, casos e controles) e também relacionar eventos e atributos (por exemplo, a densidade de eventos por população).

A razão de Kernel é a divisão da estimativa de densidade de kernel pela estimativa de densidade de kernel ponderada por um atributo. A fórmula é a seguinte:

$$\hat{\lambda}(s) = \frac{ \sum_{i=1}^{n} \frac{1}{\tau^2} k \left( \frac{s - s_i}{\tau} \right) }{ \sum_{i=1}^{n} \frac{1}{\tau^2} k \left( \frac{s - s'_j}{\tau} \right) y_i }$$

Dessa forma, temos que:

1. Numerador: Representa o alisamento dos eventos por unidade de área, considerando a função Kernel $k$ para distribuir a densidade ao redor dos pontos de interesse $s_i$;

2. Denominador: Representa o alisamento da covariável $y_i$ (por exemplo, população) por unidade de área, utilizando a mesma largura de banda $\tau$ e função Kernel $k$;

3. $\tau$: Parâmetro de suavização (largura de banda) que controla o grau de influência dos pontos vizinhos na estimativa.

4. $y_i$: Atributo associado a cada ponto $i$, como o número de habitantes ou outro fator relevante.

5. $s_i$: Localização dos eventos de interesse.

Vamos praticar no R? Nossa pequena tarefa, será criar uma “taxa suavizada”, produzir duas superfícies suavizadas e dividi-las, ponto a ponto. Para fixar o conceito, vamos montar dois conjuntos de eventos no R:

• casos: pontos gerados de forma clusterizada, simulando uma doença realmente concentrada em alguns focos;

• controles: pontos distribuídos quase ao acaso, servindo apenas de “fundo” populacional.

O resultado é a razão de Kernel, que reflete diferenças relativas entre os dois padrões:

• Valores maiores que 1 indicam maior densidade relativa no padrão agregado.

• Valores menores que 1 indicam menor densidade no padrão agregado em relação ao aleatório.

Acompanhe o código abaixo:

bandwidth <- 0.10
kernel_random <- density(aleatorio, sigma = bandwidth)
kernel_cluster <- density(agregado, sigma = bandwidth)

clu_alea_ratio <- kernel_cluster/kernel_random 

par(mfrow = c(1, 3))
plot(kernel_cluster, main = "Cluster")
plot(agregado, add = T, col = "white", pch = 20)

plot(kernel_random, main = "Aleatório")
plot(aleatorio, add = T, col = "white", pch = 20)

plot(clu_alea_ratio, main = "Razão")

Na Figura 14, note que o mapa da razão de Kernel (imagem da direita) conserva, em amarelo e vermelho, as regiões em destaque que apareciam como “áreas quentes” no kernel de casos (imagem da esquerda). Isso indica que esses agrupamentos permanecem evidentes mesmo depois da padronização. Em outras palavras, a razão de Kernel permite identificar padrões relativos e destacar áreas de concentração que seriam menos visíveis apenas com as densidades absolutas.

Agora, vamos conhecer os processos pontuais de segunda ordem, que investigam a interação entre os eventos como, por exemplo, a Função K e G. Vamos lá?

Função G - Distância do vizinho mais próximo

A função G é uma ferramenta estatística que mede a distância do vizinho mais próximo (ou seja, o evento mais próximo) em relação a um ponto específico. É útil na vigilância em saúde por ser aplicada na detecção de padrões espaciais, como agrupamentos ou dispersões.

O método do vizinho mais próximo é utilizado para estimar a função de distribuição cumulativa $\hat{G}(r)$, que se baseia nas distâncias $r$ entre eventos dentro de uma região de análise.

A função de distribuição cumulativa pode ser estimada empiricamente pela fórmula:

$$\hat{G}(r) = \dfrac{\text{#}(d(u_{i},u_{j}) \leq r)}{n}$$

Para a qual:

• $d(u_{i},u_{j})$ é a distância entre dois eventos $u_i$ e $u_j$.

• $r$ é o limite de distância considerado.

• $n$ é o número total de eventos na região de análise.

A análise gráfica de $\hat{G}(r)$ serve como uma ferramenta exploratória para identificar a existência de interações espaciais entre eventos. Os seguintes padrões podem ser observados:

• Um crescimento rápido da função para pequenas distâncias $r$ pode indicar interações entre eventos, sugerindo agrupamento ou padrões de aglomeração em escalas menores.

• Valores baixos de $\hat{G}(r)$ para distâncias pequenas, seguidos de crescimento lento, podem sugerir uma distribuição mais regular e menos interativa entre os eventos.

Os gráficos da Figura 15 apresentam as estimativas da função de distribuição $\hat{G}(r)$ para os seguintes padrões espaciais: Aleatório, Regular e Agregado. Cada gráfico compara o comportamento empírico de $\hat{G}(r)$ com a função teórica $G_{pois}(r)$, que representa um padrão de distribuição aleatório (processo de Poisson). Abaixo está a interpretação de cada caso:

1. Aleatório (Primeiro Gráfico): O comportamento de $\hat{G}(r)$ é muito próximo da função $G_{pois}(r)$, sugerindo que os eventos não apresentam interação espacial e estão distribuídos de forma aleatória. Pequenas variações são esperadas devido a flutuações empíricas.

2. Regular (Segundo Gráfico): A curva empírica $\hat{G}(r)$ cresce mais lentamente em comparação com $G_{pois}(r)$, indicando uma distribuição mais espaçada dos eventos. Esse padrão é característico de processos regulares, nos quais os eventos mantêm uma distância mínima entre si, sugerindo repulsão.

3. Agregado (Terceiro Gráfico): A curva $\hat{G}(r)$ cresce mais rapidamente para valores baixos de $r$ em relação a $G_{pois}(r)$, indicando que os eventos estão agrupados ou formam aglomerações em pequenas distâncias. Este padrão é típico de processos que apresentam interações positivas ou atração entre os eventos.

Embora o método do vizinho mais próximo seja útil para fornecer uma indicação inicial sobre a distribuição espacial dos eventos, ele é limitado a pequenas escalas espaciais, o que pode restringir a análise de padrões em distâncias maiores. Para obter uma visão mais abrangente e efetiva do padrão espacial em escalas maiores, recomenda-se o uso da função k, que permite uma análise mais detalhada ao considerar uma gama mais ampla de distâncias e a interação entre os eventos em diferentes escalas.

Função K de Ripley (ou apenas função K)

A função $K$, também conhecida como medida de momento de segunda ordem reduzida, é uma ferramenta amplamente utilizada na análise de padrões espaciais. Ela é definida para um processo univariado como:

$$\lambda K(h) = E(\text{#}\text{ de eventos contidos até uma distância } h \text{ de um evento arbitrário}),$$

Onde:

• $\lambda$ é a intensidade ou o número médio de eventos por unidade de área, assumido constante na região de estudo.

• $E()$ é o operador de expectativa.

• $h$ é a distância considerada.

A fórmula estimada para $K(h)$ é dada por:

$$\hat{K}(h) = \frac{A}{n^2} \sum_{i=1}^n \sum_{j \neq i} w_{ij}^{-1} I[x_j : d(x_i, x_j)],$$

Onde:

• $A$ é a área da região de análise.

• $n$ é o número total de eventos.

• $w_{ij}$ é o fator de correção de borda.

• $d(x_i, x_j)$ é a distância entre os pontos $x_i$ e $x_j$.

A função $K$ mensura quantos eventos estão distribuídos em círculos concêntricos ao redor de um ponto de referência (chamado de ponto focal). Esses círculos partem de um raio inicial igual a zero e aumentam gradativamente até cobrir toda a área de estudo.

Testes Genéricos de Detecção de Clusters

Os testes genéricos são aplicados quando não existe uma hipótese prévia sobre a localização de um possível cluster. Esses testes analisam todo o espaço de estudo à procura de áreas com concentração incomum de eventos. Entre os testes genéricos mais utilizados, destacamos:

i) Estatística de varredura de Kulldorff (SaTScan)

Este método utiliza janelas móveis de diferentes tamanhos para identificar regiões com alta incidência de eventos. Pode ser aplicado a modelos estatísticos diversos, como Poisson, Bernoulli ou Normal.

Exemplo prático: detectar áreas onde ocorrem surtos inesperados de doenças, sem que haja uma suposição prévia sobre onde esses clusters possam estar.

ii) Índice I de Moran

O Índice I de Moran mede a autocorrelação espacial global, verificando se há padrões significativos de agrupamento espacial no conjunto de dados. Embora não indique exatamente a localização dos clusters, ele permite confirmar se existe dependência espacial.

Exemplo prático: identificar padrões globais de altas taxas de mortalidade em uma cidade. Voltaremos a discutir esse método mais adiante no curso.

iii) Teste de Knox

O teste de Knox é especialmente útil para detectar clusters espaço-temporais. Ele avalia se eventos próximos no espaço também ocorrem próximos no tempo, indicando padrões que não são resultado do acaso, como epidemias ou surtos.

Esse teste parte da hipótese nula de que as ocorrências dos eventos são independentes tanto no espaço quanto no tempo. Caso contrário, se existir associação espacial e temporal significativa, isso sugere fortemente a presença de um cluster espaço-temporal.

Exemplo prático: identificar se casos de uma epidemia apresentam uma concentração significativa tanto no espaço quanto no tempo.

Testes Focados de Detecção de Clusters

Ao contrário dos testes genéricos, os testes focados são aplicados quando já existe um local específico de interesse e pretende-se verificar se há um cluster ao redor dele. Esse tipo de análise é indicado especialmente quando existe uma hipótese prévia sobre onde o cluster pode ocorrer. Entre os testes focados, destacamos:

i) Estatística de Varredura Focada de Kulldorff (SaTScan Focado)

Semelhante à versão genérica, esse método testa especificamente a presença de clusters em torno de pontos pré-determinados.

Exemplo prático: investigar se há aumento significativo de casos de câncer ao redor de uma usina nuclear.

Baixando e preparando os dados

library(tidyverse)
library(spatstat)
library(splancs)

# Lendo os bancos que estão no repositório github
local <- "https://raw.githubusercontent.com/ogcruz/dados_eco_2023/main/dados/"

homic <- read.table(paste0(local, "homic.dat"), col.names = c("x", "y"))
suic <- read.table(paste0(local, "suic.dat"), col.names = c("x", "y"))
acid <- read.table(paste0(local, "acid.dat"), col.names = c("x", "y"))

# Plotando os casos de homicídios em um plano cartesiano

g <- ggplot(homic) + 
  geom_point(aes(x, y, color = "Homicídios"), shape = 1) + 
  labs(color = "")
g

Porto Alegre é uma cidade disposta ao longo do eixo norte/sul. O gráfico perdeu a estrutura espacial, ajustando para o tamanho e forma da janela. Por isso é necessário informar ao programa que este tipo de objeto é um objeto espacial e tem uma escala, que deve ser preservada. Vamos fazer isso com transformando-o em um objeto espacial do pacote sf.

library(sf)

# Transformando os pontos em geometria espacial de pontos

homic_sf <- homic |> 
  st_as_sf(coords = c("x", "y"))
suic_sf <- suic |>
  st_as_sf(coords = c("x", "y"))
acid_sf <- acid |>
  st_as_sf(coords = c("x", "y"))

# Repetindo o mesmo gráfico. Repare que agora usamos a função `geom_sf()`:
g <- ggplot() + 
     geom_sf(data = homic_sf, 
             aes(geometry = geometry, color = "Homicídios"), 
             shape = 1) + 
     labs(color = "")
g

Agora vamos adicionar os pontos referentes a suicídios e acidentes:

g <- g + 
  geom_sf(data = suic_sf, 
          aes(geometry = geometry, color = "Suícidios"), 
          shape = 1) + 
  geom_sf(data = acid_sf, 
          aes(geometry = geometry, color = "Acidentes"), 
          shape = 1)
g

Agora vamos importar o polígono que corresponde ao contorno de Porto Alegre:

# Contorno de Porto Alegre
contorno.poa <- read.table(paste0(local, "/contpoa.dat"),
                           col.names = c("x", "y"))

# Plotando com a função `geom_polygon`:
g + 
  geom_polygon(data = contorno.poa, aes(x, y), 
               fill = "transparent", color = "black")

Os pontos fora do contorno são das ilhas, não devem ser incorporados a análise. Vamos transformar o contorno também em um objeto espacial (st_polygon) para identificar os pontos fora do contorno:

# Para transformar em polígono espacial,
# transformamos primeiro em matriz e em seguida
# em uma lista:

poa_sf <- contorno.poa |>
    as.matrix() |>
    list() |>
    st_polygon()

Agora sim, podemos utilizar a função st_within() para identificar os pontos fora do contorno:

# A função st_within() tem essa funcionalidade
homic_sf <- homic_sf |>
    mutate(dentro = lengths(st_within(homic_sf, poa_sf)))

ggplot(homic_sf, aes(geometry = geometry)) + 
  geom_sf(aes(color = as.factor(dentro))) +
  geom_sf(data = poa_sf, fill = "transparent") +
  labs(color = "Dentro do polĩgono")

# Filtra somente as observações dentro do polígono
homic_sf2 <- homic_sf |>
    dplyr::filter(dentro == 1)

ggplot(homic_sf2, aes(geometry = geometry)) +
  geom_sf(aes(color = as.factor(dentro))) +
  geom_sf(data = poa_sf, fill = "transparent") +
  labs(color = "Dentro do polĩgono")

Pronto, agora podemos começar a reproduzir algumas análises, que foram vistas na parte teórica a respeito dos padrões pontuais, utilizando o contexto casos de homicídios, suicídios e acidentes de carro em Porto Alegre/RS.

Primeiro, para organizar imagens e apresentá-las juntas, vamos utilizar o pacote gridExtra:

# se não estiver instalado, rodar:
install.packages("gridExtra")
library(gridExtra)

Em seguida, vamos simular alguns padrões dos dados de ponto da ocorrência de homicídios, utilizando como referência o polígono referente ao contorno de Porto Alegre/RS:

set.seed(123)

# Gerando um padrão aleatório
alea_sf <- st_sample(poa_sf, size = 100, type = "random")

g_alea <- ggplot(alea_sf, aes(geometry = geometry)) +
  geom_sf() + 
  geom_sf(data = poa_sf, fill = "transparent") +
  ggtitle("Distribuição aleatória")

# Gerando um padrão regular
uni_sf <- st_sample(poa_sf, size = 100, type = "regular") %>%
  st_as_sf() %>%
  # a função jitter é bastante utilizada e
  # serve para adicionar uma flutuação (desvio)
  # aleatória aos pontos
  st_jitter(amount = 100) %>%
  dplyr::filter(lengths(st_within(., poa_sf)) == 1)

g_uni <- ggplot(uni_sf, aes(geometry = x)) + geom_sf() +
  geom_sf(data = poa_sf, aes(geometry = geometry),
          fill = "transparent") + ggtitle("Distribuição regular")


# Gerando um padrão de cluster
# 1. Criar pontos centrais para os clusters
clusters_centers <- st_sample(poa_sf, size = 5, type = "regular") %>% 
  st_as_sf()

# 2. Adicionar pontos ao redor dos centros dos clusters
cluster_sf <- clusters_centers %>%
  st_coordinates() %>%
  as.data.frame() %>%
  rename(x = X, y = Y) %>%
  mutate(cluster_id = 1:nrow(.)) %>%
  # Adicionar 20 pontos ao redor de cada centro
  group_by(cluster_id) %>%
  do({
    center_x <- .$x
    center_y <- .$y
    n_points <- 20
    tibble(
      x = rnorm(n_points, mean = center_x, sd = 50),
      y = rnorm(n_points, mean = center_y, sd = 50)
    )
  }) %>%
  ungroup() %>%
  st_as_sf(coords = c("x", "y"), crs = st_crs(poa_sf)) %>%
  dplyr::filter(lengths(st_within(., poa_sf)) == 1) # Filtrar apenas pontos dentro do polígono

# Visualizar o padrão de cluster
g_cluster <- ggplot(cluster_sf, aes(geometry = geometry)) +
  geom_sf() +
  geom_sf(data = poa_sf, aes(geometry = geometry),
          fill = "transparent") +
  ggtitle("Distribuição em Clusters")

grid.arrange(g_alea, g_uni, g_cluster, ncol = 3)

Testando a Completa Aletoriedade Espacial (CSR)

# Convertendo para a class ppp

alea_ppp <- alea_sf |>
    as.ppp()

uni_ppp <- uni_sf |>
    as.ppp()

cluster_ppp <- cluster_sf |>
    as.ppp()

# Construindo os quadrantes com as respectivas
# contagens
alea_qc <- quadratcount(alea_ppp, nx = 5, ny = 6)
uni_qc <- quadratcount(uni_ppp, nx = 5, ny = 6)
cluster_qc <- quadratcount(cluster_ppp, nx = 5, ny = 6)

quadrat.test(alea_qc)

Chi-squared test of CSR using quadrat counts

data:
X2 = 53, df = 29, p-value = 0.008416 alternative hypothesis: two.sided

Quadrats: 5 by 6 grid of tiles

quadrat.test(uni_qc)

Chi-squared test of CSR using quadrat counts

data:
X2 = 34.789, df = 29, p-value = 0.4231 alternative hypothesis: two.sided

Quadrats: 5 by 6 grid of tiles

quadrat.test(cluster_qc)

Chi-squared test of CSR using quadrat counts

data:
X2 = 381.39, df = 29, p-value < 2.2e-16 alternative hypothesis: two.sided

Quadrats: 5 by 6 grid of tiles

Lembrando que:

• Hipótese nula ($H_0$): Os pontos estão distribuidos aleatoriamente no espaço.

• Hipótese alternativa ($H_1$): Os pontos formam agrupamentos (clusters) ou estão organizados de maneira dispersa no espaço.

Para padrão aleatório dos dados pontuais, o $p-valor = 0,008416$, ou seja, o $p-valor > \alpha$, supondo um $\alpha = 0,05$, indica que há evidências suficientes para rejeitar a hipótese nula de que o padrão de pontos segue uma distribuição aleatória completa (CSR).

Já supondo o padrão regular dos dados pontuais, o $p-valor = 0,4231$, ou seja, também o $p-valor > \alpha$, supondo um $\alpha = 0,05$, indica que não há evidências suficientes para rejeitar a hipótese nula de que o padrão de pontos segue uma distribuição aleatória completa (CSR).

Por fim, para esse cenário aglomerado dos dados pontuais, o $p-valor = 2,2 \times 10^{-16}$, ou seja, também o $p-valor < \alpha$, supondo um $\alpha = 0,05$, indica que há evidências fortes para rejeitar a hipótese nula de que o padrão de pontos segue uma distribuição aleatória completa (CSR).

Processo de primeira ordem: Gerando a estimativa de densidade de Kernel (mapa de calor)

Agora iremos verificar o padrão espacial de primeira ordem dos casos dos casos de homicídios, suicídios e acidentes de carro em Porto Alegre/RS.

homic_ppp <- homic_sf2 |>
  as.ppp()
Window(homic_ppp) <- as.owin(poa_sf)

homic_den1 <- density(homic_ppp, sigma = 200, diggle = TRUE)

g_homic_den1 <- ggplot(as_tibble(homic_den1), aes(x, y)) + 
  geom_tile(aes(fill = value)) +
  geom_point(data = as_tibble(homic_ppp), 
             color = "white", shape = 1, size = 0.5) + 
  geom_sf(data = poa_sf, aes(geometry = geometry), 
          fill = "transparent", inherit.aes = FALSE) + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Homicídios", subtitle = "sigma = 200") +
  guides(fill = "none") + 
  theme_void()

suic_ppp <- suic_sf |>
  as.ppp()
Window(suic_ppp) <- as.owin(poa_sf)

suic_den1 <- density(suic_ppp, sigma = 200, diggle = TRUE)

g_suic_den1 <- ggplot(as_tibble(suic_den1), aes(x,y)) + 
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(suic_ppp),
             color = "white", shape = 1, size = 0.5) + 
  geom_sf(data = poa_sf, aes(geometry = geometry), 
          fill = "transparent", inherit.aes = FALSE) + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Suicídios", subtitle = "sigma = 200") +
  guides(fill = "none") + 
  theme_void()

acid_ppp <- acid_sf |>
  as.ppp()
Window(acid_ppp) <- as.owin(poa_sf)

acid_den1 <- density(acid_ppp, sigma = 200, diggle = TRUE)

g_acid_den1 <- ggplot(as_tibble(acid_den1), aes(x, y)) + 
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(acid_ppp),
             color = "white", shape = 1, size = 0.5) + 
  geom_sf(data = poa_sf, aes(geometry = geometry),
          fill = "transparent", inherit.aes = FALSE) + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Acidentes", subtitle = "sigma = 200") +
  guides(fill = "none") + 
  theme_void()

grid.arrange(g_homic_den1, g_suic_den1, g_acid_den1, ncol = 3)

Agora iremos investigar verificar o padrão espacial de segunda ordem dos casos dos casos de homicídios em Porto Alegre/RS utilizando diferentes larguras de banda ($\sigma$ = 100, 200 e 500 metros).

homic_den2 <- density(homic_ppp, sigma = 100, diggle = TRUE)

g_homic_den2 <- ggplot(as_tibble(homic_den2), aes(x, y)) + 
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(homic_ppp),
             color = "white", shape = 1, size = 0.5) +
  geom_sf(data = poa_sf, aes(geometry = geometry), 
          fill = "transparent", inherit.aes = FALSE) + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Homicídios", subtitle = "sigma = 100") +
  guides(fill = "none") + 
  theme_void()

homic_den3 <- density(homic_ppp, sigma = 500, diggle = TRUE)

g_homic_den3 <- ggplot(as_tibble(homic_den3), aes(x, y)) + 
  geom_tile(aes(fill = value)) + 
  geom_point(data = as_tibble(homic_ppp),
             color = "white", shape = 1, size = 0.5) + 
  geom_sf(data = poa_sf, aes(geometry = geometry), 
          fill = "transparent", inherit.aes = FALSE) + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Homicídios", subtitle = "sigma = 500") +
  guides(fill = "none") + 
  theme_void()

grid.arrange(g_homic_den2, g_homic_den1, g_homic_den3, ncol = 3)

Vamos fazer a razão de kernel entre as causas de homicídio e suicídios de Porto alegre/RS:

suic_homic_ratio <- suic_den1/homic_den1

g_suic_homic <- ggplot(as_tibble(suic_homic_ratio), aes(x, y)) + 
  geom_tile(aes(fill = value)) + 
  geom_sf(data = poa_sf, aes(geometry = geometry), 
          fill = "transparent", inherit.aes = FALSE) + 
  scale_fill_viridis_c(option = "B") +
  labs(title = "Razão", subtitle = "Suicídios/Homicídios") +
  guides(fill = "none") + 
  theme_void()

grid.arrange(g_homic_den1, g_suic_den1, g_suic_homic, ncol = 3)

Processo de segunda ordem: Funções G e K

Verificando o padrão espacial de segunda ordem dos casos de homicíodos em POA/RS:

# se não estiver instalado, rodar:
install.packages("splancs")
library(splancs)

library(splancs)

homic_ppp2 <- as.ppp(homic_sf2$geometry)

plot(homic_ppp2, pch = 19, cex = 0.5)
polymap(contorno.poa, add = T)

par(mfrow = c(1, 2))
plot(envelope(Y = homic_ppp2, fun = Gest, nsim = 9), main = "Funcao G")

Generating 9 simulations of CSR … 1, 2, 3, 4, 5, 6, 7, 8, 9.

Done.

plot(envelope(Y = homic_ppp2, fun = Kest, nsim = 9), main = "Funcao K")

Generating 9 simulations of CSR … 1, 2, 3, 4, 5, 6, 7, 8, 9.

Done.

Os dois gráficos (funções G e K) sugerem que os homicídios em Porto Alegre (RS) apresentam um padrão agrupado em vez de uma distribuição aleatória. Esse padrão indica que há concentrações de homicídios em determinadas áreas, o que pode estar associado a fatores socioeconômicos, geográficos ou outros. Para uma análise mais aprofundada, seria necessário investigar as causas do agrupamento.